Importar y configurar repositorio evolution_of_species

Estas celdas es para importar la librería desde GitHub y que se mantenga actualizado. Hay que descomentar las líneas si se utiliza desde google colab.

Librerías

Constantes

Funciones que encapsulan integrales

Se definen las funciones \begin{equation} F(t, x) = \int K(x, y) R(t, y) dy \end{equation} y \begin{equation} G(t, y) = \int r(x) K(x, y) u(t, x) dx \end{equation}

El proposito de las clases IntegrarF e IntegrarG es proporcionar una forma sencilla de realizar estos cálculos de forma eficiente.

Para esto representa estas funciones por medio de matrices, asumiendo que $K$, $R$, $r$ y $u$ vienen en forma matricial. Además de que una vez fijadas $K$ o $r$, luego no será necesario llamarlas en cada instante; sólo será necesario actualizar $R$ y $u$.

La forma de utilizar es definiendo instancias de la forma

F = FunctionarF(K, R_0, x_lims=(0, 1), y_lims=(0, 1), dt=0.01, N=100, M=100, T=100)
G = FunctionarG(r, K, u_0, **PARAMETROS_DISCRETIZACION)

donde los parámetros de discretización son opcionales

Para calcular los valores basta con preguntar la componente $(n, j)$-ésima y lo cálcula al mismo tiempo que lo entrega. También se pueden actualizar las matrices $R$ y $u$ en una misma línea. Además, se puede consultar por la matriz que lo representa. Por ejemplo:

# Retorna el valor de una componente
F[3, 4]
# Retorna el valor de una fila
G[3]
# Actualiza la matriz R y retorna la fila
F.actualize_row(R[3], 3)[3]
# Consulta por su matriz
G.matrix

Solvers para calcular R

Dado que tenemos dos métodos para calcular las Ecs. de Perthame, en donde la única diferencia radica en la actualización de $R$, se diseñaran dos solver que encapsulen las distintas formas de calcularlo.

El primer método consiste en resolve la EDO \begin{equation} \partial_t R(t, y) = -m_2(y) R(t, y) + R_{in}(y) - R(t, y) G(t, y) \end{equation} con $G(t, y)$ de la forma que se comentó anteriormente.

El segundo método consiste en asignar a $R$ el siguiente valor: \begin{equation} R(t, y) = \frac{R_{in}(y)}{m_2(y) + G(t, y)} \end{equation}

El primer solver utiliza el siguiente esquema de actualización \begin{equation} R^{n+1} = (1 - \Delta t \cdot (m_2 + G^{n})) \odot R^{n} + \Delta t R_{in} \end{equation}

donde $\odot$ representa la multiplicación componente a componente.

El segundo solver utiliza el siguiente esquema de actualización \begin{equation} R^{n+1} = \frac{ R_{in} }{ m_2 + G^{n+1} } \end{equation}

Y para el tercer solver se utilizó una aproximación centrada de segundo orden para la derivada temporal, resultando en el siguiente esquema: \begin{equation} R^{n+1} = R^{n-1} - 2 \Delta t \cdot (m_2 + G^{n}) \odot R^{n} + 2 \Delta t R_{in} \end{equation}

Para ocuparlo basta con generar una instancia de alguno de los solvers:

solver = Solver1R(r, K, u_0, R_0)

donde los parámetros de actualización son opcionales.

Al igual que antes, se puede actualizar u. Además, se puede buscar la componente $(n, k)$-ésima de la matriz.

# Calcula y retorna la fila n
solver[n]
# Calcula y retorna la componente n, k
solver[n, k]
# Actualiza u y retorna la componente
solver.actualize_row(u[2], 2)[n, k]
# Actualiza el solver al paso n + 1
solver.actualize_step_np1(1)
# Calcula la masa total (en este caso, los recursos totales) a lo largo del tiempo
solver.calculate_total_mass()
# Retorna una función R(t, y) interpolada
R = solver.function_interpolated()

Solvers para calcular $u(t,x)$

La ecuación

$$ \partial_t u(t,x) = u(t,x)(-m_1(x) + r(x)F(t,x)) + \epsilon \Delta u(t,x) $$

se puede reescribir usando $A(t,x) = -m_1(x) + r(x)F(t,x)$ y se ve de la siguiente forma:

$$ \partial_t u(t,x) - \epsilon \Delta u(t,x) = A(t,x) u(t,x) $$

Que es muy similar a la ecuación del calor, salvo que en este caso $f(t,x) = A(t,x)u(t,x)$ depende de $u$, y el Laplaciano está siendo perturbado por un término $\epsilon$. Es por ello, que se adaptará el esquema de discretización de la ecuación del calor para resolver este problema, incorporando la dependencia de $u$ de la función $f$ y la pertubación.

Entonces, reemplazando en la ecuación con las diferencias finitas forward de primer y centrada de segundo orden, se obtiene:

$$ \frac{u^{n+1}_j - u^{n}_j}{\Delta t} - \epsilon \frac{u^{n}_{j-1} - 2u^{n}_j - u^{n}_{j-1}}{h^2} = A^n_j u^n_j $$

que al despejar $u_j^{n+1}$ es de la siguiente forma:

$$ u_j^{n+1} = \alpha \ u_{j-1}^{n} + \beta_{j}^{n} \ u_{j}^{n} + \gamma \ u_{j+1}^{n} $$

con $\alpha$, $\beta$ y $\gamma$ $$ \alpha = \epsilon \frac{\Delta t}{h^2} \qquad \beta_{j}^{n} = 1 - 2\alpha + \Delta t A^n_j \qquad \gamma = \alpha $$

Resolver ecuación Perthame

Teorema de condición de supervivencia

Con las hipótesis necesarias, si $\int \frac{R_{in}(y)}{m_2(y)}|\ln R^0(y)|dy < \infty$ y se tiene (para el caso Gaussiano) \begin{equation} \bar{m}_1 \bar{m}_2 \geq \frac{ r M_{in} }{ \sqrt{2\pi(\sigma_K^2 + \sigma_{in}^2)} } \end{equation} Entonces la solución se extinguirá, es decir, que $\int u(t, x) dx$ converge a cero mientras que $R(t, y)$ converge.

Parámetros para la extinción

Se utilizaron $\varepsilon=0.001$, $\sigma_K=0.6$, $\sigma_{in}=1$, $r=1$, $M_{in}=3$ y $\bar{m}_1=\bar{m}_2=1.1$

Caso cuasi-estático

Parámetros para la supervivencia

Se utilizaron $\varepsilon=0.001$, $\sigma_K=0.6$, $\sigma_{in}=1$, $r=1$, $M_{in}=3$ y $\bar{m}_1=\bar{m}_2=1$

Caso cuasi-estático

Divergencia evolutiva

Se eligieron:

Caso cuasi-estático

Otros casos